Dynamic Kosterlitz-Thouless transition in 2D Bose mixtures of ultra-cold atoms 
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We propose a realistic experiment to demonstrate a dynamic Kosterlitz-Thouless transition in 
ultra-cold atomic gases in two dimensions. With a numerical implementation of the Truncated 
Wigner Approximation we simulate the time evolution of several correlation functions, which can 
be measured via matter wave interference. We demonstrate that the relaxational dynamics is well- 
described by a real-time renormalization group approach, and argue that these experiments can 
guide the development of a theoretical framework for the understanding of critical dynamics. 
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The understanding of non-equilibrium phenomena, in 
particular dynamic phase transitions, is an open fron- 
tier in many-body physics. While for equilibrium sys- 
tems a rich variety of methods has been established, non- 
equilibrium systems are notoriously difficult to grasp, and 
a functioning conceptual framework is lacking. Given this 
state of research, ultra-cold atomic systems can play a 
crucial role in further understanding many-body dynam- 
ics. In fact, the unprecedented control of well-defined, 
isolated systems of ultra-cold atoms, has led them to 
be considered as 'quantum simulators' [l|: Cold atom 
systems are manipulated to create paradigmatic model 
systems of condensed matter physics, such as the Bose- 
Hubbard model [2.], spin chains [3], the unitary Fermi 
gas magnetic systems [5[, the Dirac equation j6[ and 
equilibration in one-dimensional (ID) gases f^. The ex- 
perimental measurements are cross-checked with theory, 
to guide the development of a theoretical framework [8[. 

In this paper, we propose to apply the concept of quan- 
tum simulation for a non-equilibrium setup, in particular 
for dynamic phase transitions. We present an implemen- 
tation of the time-dependent renormalization group (RG) 
description, derived in [9], and demonstrate that it quan- 
titatively predicts the dynamic evolution across a critical 
point, by comparing it to a numerical simulation. In par- 
ticular, this proves the universality of the relaxational 
dynamics in this system, a concept well-established for 
equilibrium phase transitions, however not developed for 
dynamics in closed systems. Other RG treatments of 
non-equilibrium systems were reported in [lo]. 

We consider a system of weakly interacting bosons in 
two dimensions (2D), for a review see e.g. [11], which in 
equilibrium undergoes a Kosterlitz-Thouless (KT) phase 
transition as a function of temperature . This system 
has two thermal phases, defined through the long-range 
behavior of the two-point correlation function G(r) = 
(?/^^(0)?/^(r)), where ?/^(r) is the particle annihilation op- 
erator at site r. At low temperatures this function decays 
algebraically, G(r) ~ |r|~'^/^. The exponent r increases 
monotonically from zero to 1 as the temperature is in- 
creased from zero to the critical temperature Tc. Above 
Tc, the functional form of G(v) changes to exponential 




FIG. 1: We prepare a 2D atom cloud in state 1 (blue), and 
apply a 7r/2 pulse at We apply a field gradient at ts, 

which separates state 1 and 2 (red) spatially. We release the 
atoms at time tr and measure their interference properties. 

decay, G(r) ~ exp(— |r|/ro), with some decay length 
tq. This change is due to the deconfinement of vortex- 
antivortex pairs, and defines the KT transition. We trig- 
ger this transition dynamically by a quench, as described 
below. We find that after an intermediate time the sys- 
tem develops a metastable state in which the phononic 
modes have equilibrated, and G(r) shows algebraic scal- 
ing with an exponent r, which can be larger than the 
critical value 1. This exponent then increases slowly in 
time, until the correlation function changes to exponen- 
tial scaling, indicating dynamic vortex deconfinement. 
We demonstrate that the real-time evolution of r(t) can 
be described by a real-time RG approach, by comparing 
it to our simulations. This constitutes a conceptually new 
insight into many-body dynamics. For the ultimate vali- 
dation of this theoretical approach we propose a specific 
experimental setup in this paper. 

In Fig. [1] we sketch the quench and measurement se- 
quence. We consider a 2D gas of atoms, such as ^^Rb, 
at a temperature below Tc with an initial scaling expo- 
nent < 1 in an internal state labelled 1. At some 
time ^7^/2, a radio- frequency or microwave pulse drives 
the transition from state 1 to another state 2. The dura- 
tion and intensity of the pulse are adjusted to provide a 
7r/2 transition, so that the superfluid (SF) is now in the 
state ('0i(r) + ?/;2(r))/>/2, where V^i(r) and ^^2(1*) are the 
single-particle operators of states 1 and 2. If the inter- 
action strengths between 1 and 2 are identical (as they 
approximately are for the hyperfine levels of ^^Rb in its 
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FIG. 2: a) gi{x,t)/no on a logarithmic scale, for U/J = 0.25 and n = 0.89 on a 250 x 250 lattice. The 7v/2 pulse is applied 
at Jt7r/2 = 1000. Ui2 is turned off at Jts = 1200. After ts, gi(x,t) gradually changes its functional form, from algebraic to 
exponential, b) Fitted algebraic exponent r obtained from the data in a), c) Fitting errors Ee and Ea- 



ground state), the rotation between the internal states is 
a symmetry of the system and the resulting SF state is 
still a steady state. Then at a later time tg, we perform a 
quench of the gas by spatially separating the two species 
with a magnetic field gradient. This has the effect that 
the inter-species interaction is set to zero. For ^^Rb, the 
two hyperfine states could be \F = I^tuf = —1) and 
|F = 2,mF = — 1), which have Lande factors of opposite 
sign. We measure the coherence properties at some re- 
lease time tr by (i) switching of? the confinement along 
the initially frozen direction, so that the two gases expand 
along this direction and overlap; (ii) performing another 
7r/2 pulse between 1 and 2; (iii) recording the atomic 
density in one of the two states, say 1, which reveals the 
matter-wave interference between the two planes. We 
note that 'splitting' of a ID SF was discussed in [13]. 

We describe the system with the Hamiltonian H = 
i^i + + Hi2, where 



(1) 



and 



Hi2 = [512^1 (r)^|(r)^2(r)^i(r)], (2) 



where the states are labeled a = 1,2, and the dimension 
D = 2] I is the discretization length scale of this repre- 
sentation, see Ref. [14], m is the atom mass, go {912) is 
the strength of the intra- (inter-) species contact inter- 
action. For atoms confined to 2D motion by a harmonic 
potential in the third direction, these are approximately 
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lo is the harmonic oscillator length lo = {Ti/muJoY^'^ of 
the confining potential muJoZ^ /2 in z direction; and 
ai2 are the s-wave scattering lengths. For hyperfine states 
of ^^Rb these are around 5 nm. 



Initially all atoms are in state a = 1, and form a SF 
of total density po at a temperature T < Tc. We de- 
scribe this state with a Bogoliubov analysis b ased on the 
phase-density representation ^^(r) ^ g*6'(r)^^^ _^ Sp{r), 
see Ref. [14]. As mentioned, the correlation function 
Gi{r) = {ipl{0)ipi{r)) decays algebraically, G'i(r) ~ 
\r\-^^^, with 



7rh^pJ{2mT), 



(4) 



see Refs. [H, [I^, implicitly defining the SF density p^. 
Within the Bogoliubov approximation, ps = po- How- 
ever, based on Ref. [17], for weak interactions and away 
from the critical regime, an improved estimate is 



7rh^Po/{2mT)^Co, 



(5) 



see [3, with Co = {\n{2gm/h^) - l)/4. We simu- 
late the dynamics with a numerical implementation of 
the Truncated Wigner approximation (TWA), see e.g. 
Refs. [1^, [20|. It is convenient to represent the system 
as a Hubbard model with two species where i is 
the lattice site, and a = 1,2 the species index. The 
single-particle operators of the two representations are 
related as ba^i = /^/^V^a(ri), where is the real-space 
location corresponding to lattice site i. The tunneling 
energy J and m are related through J = / {2mP)] go-, 
gi2 and the Hubbard interactions through U = go/l^ 
and U12 = 912/1^ ^ respectively. We use h/J as the time 
unit and choose U/J = U12/J = 0.25. Using the rela- 
tion U/J = V327ras//o from Eq. (|3]), this value corre- 
sponds to lo ~ 200nm, which is routinely achieved ex- 
perimentally. We sample the initial state according to 
the Wigner distribution VK(xk,Pk) ^ ^^P(~^k/^^^ k ~ 
Pk/2^p,k). with a- 



2 ^ = l/(2a;ktanh(a;k/(2To))) and 
co'k/(2 tanh(cJk/(2To))), where Tq is the initial- 
iz ation temperature, and the Bogoliubov modes /3k = 
Y^Co'k/2xk + i/^/2uJ^p\^. From these we construct Sp{r) 
and 6{r) and from these '^i(r), see |21|]. '^2(1*2) is initial- 
ized with vacuum fluctuations (|V^2(ri)P) = 1/2/^. The 
initial state is propagated by the equations of motion. 

In Fig. [2] a) we show ga{x,t) = {h\^^{t)ha^i+x{t)) for 
a = 1 for a 250 x 250 lattice, for no = l^ po = 5. We 
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FIG. 3: a) Algebraic exponent r(t) from fitting 5'i(x, t) on a 100 x 100 lattice, for different r^. The red, dotted line is the critical 
value 1. The blue, dashed line is r(t) for n = 0.89 (as in Fig. [2]). b) Algebraic and c) exponential fitting errors Ea,e- 



use To/J = 7 in the initialization, to let the system equi- 
librate to a SF with scaling exponent = 0.89. We 
apply the 7r/2 pulse at time = 1000; the den- 

sity drops to half its initial value. At time Jts = 1200 
we turn off [/12, see [23j and [23]. The evolution af- 
ter the quench separates into two time domains. The 
system first relaxes to a metastable SF state, and then 
changes the functional form of the correlation function 
from algebraic to exponential on a much longer time 
scale, t ^ lO^h/J. To study this quantitatively, we fit the 
correlation function with two functions. We use fa{x) = 
C{s'm{7rx/N)N/7T)~^^^ with C and r as fitting parame- 
ters to test for algebraic scaling. N is the number of lat- 
tice sites in one dimension. For large fa{x) ~ \x\~^^'^. 
The fitted exponent r is plotted in Fig. [2]b), and for a 
shorter time evolution and a 100 x 100 lattice as the blue 
dotted line in Fig. [3] a), see [25|. The constant value of r 
before is of our physical initial ensemble. For expo- 
nential scaling, we use fe{x) = C ex.-p{—f<ism{7Tx/N)N/7T) 
with C and hi as fitting parameters. We define the fitting 
errors Ea^e = Y.fji{9{x,t) - fa,ef- In Fig. [2] c), we 
show the fitting errors for the data plotted in Fig. [2] a). 
We recover the time scales discussed above, (i) After the 
random draw of Xk and pk and before the 7r/2 pulse, the 
system rapidly equilibrates in few tens ofh/J and the two 
errors fluctuate around a ratio of E^/Ea ^ 100, indicat- 
ing that gi (x, t) is much better fitted with an algebraic 
test function, (ii) After the 7r/2 pulse and before the 
quench, the density of species 1 drops to half its original 
value after ^7^/2 but the gas is still in equilibrium. Both 
errors Ea^e thus drop by a factor of 4, while maintaining 
the ratio ~ 100. (iii) After the quench, E^ decreases and 
Ea increases on a long time scale. 

We repeat this calculation for different r^, and plot 
rit) and Ea^e as a function of time and of r^, see Fig. [3l 
For small r^, r(t) increases initially and then remains 
constant below 1, and Ee/Ea ^ 10^. For n > 0.6, r(t) 
increases above 1, and does not remain constant, but 
shows a slow increase. Fig. [3] a). Ea grows, and becomes 
eventually larger than E^^ as shown in Fig. [3]b) and c), 
which indicates the functional form change of ^i(x, t). 
This demonstrates a dynamic KT transition. We also 



observe a time regime in which r > 1, but the change 
to exponential scaling has not yet occurred. This is the 
metastable, supercritical SF described in Refs. [il, [l5| . 
Further we note that both at initial time and after the 
second quench, the algebraic error briefly spikes up (see 
also Fig. 2 of Supp. Mat.). This is due to the light cone 
dynamics [9|, [HI, |22| triggered by these two quenches: 
After each quench, gi{x,t) is only piece- wise algebraic, 
i.e. it behaves as Ixl"'^^/^ for \x\ < c(t — tg)^ and as 
|^|-r2/4 ^ with, differing exponents ri,r2; 

thus the fitting error is increased. 

We now demonstrate that this transition can be cap- 
tured with the RG approach derived in Ref. [9|]. Its 
key statement is that the transition from the metastable, 
'pre-thermalized' state to the fully-thermalized, vortex- 
deconfined state is described by the flow equations 

1 = 2(1- l/r)y; J = M^'ay^r (6) 

where a is a non- universal constant, and y is the vor- 
tex fugacity, see and 0. I is the flow parameter, 
related to real time by t = toe^, with some time con- 
stant to- We thus have to determine the scaling expo- 
nent of the met a-st able state, referred to as r*. The 
system has two phononic sectors, corresponding to sym- 
metric and anti-symmetric superpositions of phase and 
density excitations. Immediately after the quench, only 
the symmetric sector is thermally activated, whereas 
the anti-symmetric sector only displays vacuum fluctu- 
ations. We find that the two sectors thermally equi- 
librate on a short time scale of ^ lO^fi/J, by calcu- 
lating ratios of ^s(x) = (?/^i(O)^?/^2(O)^'02(x)V^i(x)) and 
ga{-s) = {M^yM^)4i^)M^)). such as ^,(r)/^(r)2, 
^a(r)/^(r)2, gs{r)ga{r)/g{T)^, etc, which ah approach 
unity. This result differs from the behavior of ID gases, 
for which the two sectors stay out-of-equilibrium for a 
much longer time [l3|. To determine r*, we use energy 
conservation. For T much larger than the mean-field en- 
ergy, the total energy scales as T^. After the quench and 
equilibration, the resulting temperature is thus T/^/2. 
We use Eq. (|5j) for r, note that the density is reduced 
by 1/2, and find r*^ = (l/(V^rO + (1 - 1/V^) * Co)-\ 
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FIG. 4: a) Numerical data for r, averaged over Jt — [2100, 3000]. The red, filled dots have r < 1, the dark circles have r > 1. 
The black, continuous line is ti with the optimal fitting parameters, fitted to the data r < 1. The red, dotted line indicates the 
critical line r = 1. The black, dashed line is Ti*^. b)-d) Numerical data r at times ti = 1500/ J, t2 = 2400/ J and ta = 2970/ J 
(b-d) as a function of n. The blue lines labelled I - III are in b) t(^* = 1.7), t(^* = 2.1) and r(^* = 2.5), respectively. In c) 
T{t = 2.1), T{t = 2.5) and T{t = 2.9), respectively, and in d) t(/* = 2.3), T{t = 2.7) and T{t = 3.1), respectively. 



with Co = {\n{2gm/n^) - l)/4, for a 2D gas in contin- 
uum. However, Co will be renormalized for a discretized 
representation, see Ref. [17], and we use it as a fitting 
parameter. We show r^*^ in Fig. |4] in comparison to the 
numerical data, with the optimal Co determined below. 
We use Tj*^ as the initial value for r(0) = r*:^^. To deter- 
mine the initial value for the fugacity, we write the flow 
equation as d{ay'^)/dl = 4(1 — l/r){ay'^). The quantity 
(r — 1)^ — 327r^a?/^ is invariant under the flow, and thus 
the asymptotic value for r* below the critical point is 
r*(oo) = 1 - ((r - 1)2 - 327r2a7/2)i/2. This motivates to 
use Ti = 1 - ((tj*^ - 1)2 - 327r2A)^/2 ^ fitting function 
for Co and A = ay'^{0). In Fig. |l]a) we show r^*, with 
Co = 0.299 and A = 2.7 x 10~^, and the numerical data 
for r averaged over the time range Jt = [2100, 3000]. We 
use these values in our initial conditions. We integrate 
the flow equations to different values £*, in particular to 
0, 0.1, 0.2, 4.0. In Figs, m (b-d) we show the numeri- 
cal results for r at the times ti = 1500/ J, t2 = 2400/ J 
and ts = 2970/ J, each averaged over a time interval of 
60/ J. We determine the value of for which r(£*) fits 
the numerical result the closest. The ratios ti/tj are ap- 



proximately ex.-p{£i/£j), but we note that because of the 
logarithmic dependence of ^ on t, a large numerical un- 
certainty is present. We show the optimal t(^*), and two 
close-by solutions for visual comparison. We find that 
the RG flow well describes the critical dynamics. 

In conclusion, we have presented a realistic experiment 
to investigate the dynamic KT transition in ultra-cold 
gases in 2D. We demonstrate that the critical dynam- 
ics can be described by the RG approach developed in 
Ref. [9], as it predicts correctly the dynamic evolution 
found in the TWA simulation. The time evolution of 
the correlation functions can be detected via interference 
measurements discussed in Ref. pJJ]. These predictions 
and their experimental verification would pave the way 
to an RG-based theory framework for critical dynamics. 
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Appendix A: Supplemental Material 

We show the single-particle correlation function g{x^t) 
for a 100 X 100 lattice, for = 0.89 in Fig. [5l for a time 
interval of Jt = [0,3000]. The 7r/2 pulse is applied at 
-^^7r/2 = 1000, where the density drops to half, and the 
quench is applied at Jts = 1200. As described in the 
main text, we fit this correlation function with an alge- 
braic and an exponential fitting function. The algebraic 
fit gives r{t) which is depicted as a blue, dashed line in 
Fig. 3 a) of the main text. In Fig. [6] we show the two 
fitting errors for the data shown in Fig. O We see the 
same behavior as for the example in Fig. 2 in the main 
text, but due to the shorter integration time the interme- 
diate time behavior is better resolved, and because of the 
smaller system size the numerical uncertainty is reduced. 
Before the 7r/2 pulse, the ratio Ee/Ea ^ 10^, after the 
pulse each error drops to a quarter of its value while main- 
taining the same ratio. After the quench at Jts = 1200, 
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FIG. 5: Single-particle correlation function g{x,t) for n 
0.89, for a 100 x 100 lattice. 
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FIG. 6: Algebraic and exponential fitting errors Ea,e for 
0.89, for a 100 x 100 lattice. 



the ratio is initially of that same magnitude, but then 
changes on a much longer time scale. On an interme- 
diate time scale, there is a met ast able regime in which 
the correlation function is better fitted algebraically, but 
the exponent r is much larger than 1, i.e. the system is 
supercritical. 

To construct the data shown in Fig. 3 of the main text, 
we repeat calculations as shown in Fig. [5l for a 100 x 100 
lattice, and for different r^, and we determine r(t), Ea{t) 
and Ee{t) from it. 



